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ABSTRACT 


A  derivation  is  presented  of  the  basic  equations  of  the  Quasi- 
Wundy  Code;  this  code  is  widely  used  for  computer  calculations 
relating  to  explosive  ^sterns.  The  paper  begins  with  the  basic 
notions  and  attempts  to  clear  up  the  confusion  concerning 
Lagrangian  and  Eulerian  equations  which  is  all  too  common  in 
the  literature.  A  brief  derivation  of  the  fundamental  equations 
of  continuum  mechanics  follows  and  these  equations  are  then  differ¬ 
enced  to  give  the  equations  of  the  code  in  the  form  actually 
employed  in  the  code.  The  computational  procedure  is  then  discussed 
and  two  questions  which  often  anise  with  respect  to  the  wording  of 
the  code  are  answered.  It  is  hoped  that  this  paper  will  provide 
a  background  for  answering  other  questions  as  they  arise. 
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FOREWORD 


The  wide  use  of  the  Quasi  Wundy  Code  for  various  calculations 
rela  ins  to  explosive  systems  has  shown  the  need  for  an  adequate 
summary  of  the  basic  equations  and  computational  schemes  of  the 
code  together  with  the  derivations  behind  them.  To  the  best  of 
the  author's  knowledge,  such  information  has  not  been  documented 
and  is  not  readily  available  from  any  source,  although  it  is 
necessary  in  order  to  answer  some  of  the  questions  which  arise 
with  respect  to  use  of  the  code.  Answers  to  two  such  questions 
have  been  included  in  the  body  of  the  report  and  it  is  hoped, 
that  sufficient  insight  into  the  code  is  generated  to  prepare 
the  reader  to  answer  any  other  questions  which  may  arise. 

This  paper  should  serve,  further,  both  to  introduce  the 
reader  to  the  code  as  well  as  to  make  known  its  full  possibilities 
so  that  individual  investigators  can  adapt  it  to  their  own  special 
problems.  In  accordance  with  this  goal,  the  paper  proceeds  from 
basic  notions  of  fluid  mechanics  to  the  Fortran  statements  of 
the  code.  Some  material  of  a  rather  basic  nature  has  been 
included,  but  it  is  felt  that  this  is  worthwhile  since  it  leaves 
the  paper  more  nearly  self  contained  thus  ensuring  uniformity  of 
notation  and  circumventing  the  difficulty  of  obtaining  good 
reference  material. 
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I.  BASIC  NOTIONS  AND  COORDINATE  SYSTEMS 


The  Quasi  Wundy  Code  is  a  one  dimensional  computer  code 
designed  to  calculate  various  quantities  of  interest  in  an 
explosive  system.  It  is  assumed  that  all  materials  involved 
obey  certain  basic  equations  of  fluid  mechanics  which  express 
the  basic  conservation  laws:  energy,  momentum  and  mass. 

In  order  to  obtain  these  equations,  it  is  necessary  to  make 
a  few  preliminary  remarks.  The  first  concerns  the  description 
of  the  system  which  is  reacting.  There  are  two  possible  ways  of 
specifying  physical  quantities  such  as  pressure  and  density. 

The  first  is  simply  to  set  up  a  coordinate  system  and  locate  e 
point  (x,y,z)  in  that  coordinate  system.  Then  the  pressure 
(for  example)  is  defined  as  some  function  of  the  coordinates  of 
the  point,  (x,y,z)  and  the  time.  This  is  the  Eulerian  approach. 

A  second  alternative  is  to  identify  each  particle  of  fluid, 
not  by  its  present  location  relative  to  a  set  of  fixed  coordinates, 
but  by  the  location  at  time  t  =  0.  This  identifies  each  particle 
of  the  fluid  in  terms  of  its  original  location.  This  is  the 
Lagrangian  approach. 

In  both  cases,  a  set  of  fixed  axes  are  chosen  and  the 
motion  is  described  relative  to  these  fixed  axes.  The  differ¬ 
ence  lies  in  the  choice  of  independent  variables.  Euler 
chooses  the  location  in  space  and  time  and  seeks  to  determine 
pressure,  density,  etc.  at  that  point.  Lagrange  chooses  the 
initial  coordinates  of  the  "particle"  of  fluid  and  seeks  to 
determine  its  present  location,  the  pressure,  density  and  so  on. 

It  is  to  be  emphasized  that,  mathematically,  the  difference 
between  these  two  approaches  is  a  difference  in  viewpoint. 

The  solution  to  the  Lagrange  equations  include  x  =  f(a,b,c,t) 
as  the  x  coordinate  of  the  present  position  of  the  particle. 

This  is  the  Euler  coordinate,  x.  Conversely,  the  solutions  to 
the  Eulerian  equations  will  include  a  function  which  will 
describe  the  past  and  future  locations  of  the  particle  now  at 
the  point  (x,y, z).  Constants  of  integration  will  occur  which 
depend  upon  the  position  of  the  particle  at  t  =0.  These 
coordinates  are  the  Lagrange  coordinates. 

The  Lagrange  approach  is  better  suited  than  the  Euler  for 
describing  the  motion  of  mixtures  of  fluids  with  different 
densities  and  equations  of  state,  since  it  is  only  necessary  to 
assign  these  properties  to  the  initial  state  of  the  fluid;  during 
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the  subsequent  motion  it  is  quite  easy  to  determine  which  fluid 
is  described  by  some  calculated  value  of  pressure  for  example,  j 

since  the  fluid  is  still  labelled  with  its  original  coordinates 
in  the  Lagrangian  scheme. 

This  point  is  worth  developing  in  detail  since  the 
equations  of  the  code  are  developed  in  Lagrangian  form  while 
.the  ordinary  approach  to  fluid  mechanics  is  by  the  Euler 
equations.  Accordingly  we  give  an  example  to  illustrate  this 
point. 


FIGURE  1 
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In  Figure  l,x  and  y  are  the  Eulerian  coordinates  and  a  and 
b  are  the  Lagrangian  coordinates.  In  the  Lagrangian  scheme, 
the  particle  is  characterized  by  the  coordinates  a,  b  and 
retains  these  labels  for  its  entire  motion.  Its  location  at 
any  time  is  given  by  x  =  f(a,b,t)  and  y  =  g(a,b,t).  Any 
physical  quantity  such  as  density  is  determined  as  a  function 
of  a,  b,  and  t:  p  =  p(a,b,t).  This  is  the  present  density  of 
the  particle  ■which  was  originally  at  the  point  ( a,b) .  If 
Euler's  equations  are  solved,  the  solution  will  include  a 
density  equation  p  =  p/(x,y,t).  This  is  the  density  of  the 
fluid  now  at  the  point  (x,y).  if  (x,y)  is  the  present  location 
of  the  particle  originally  at  (a,b),  we  can  substitute  x  =  f(a,b,t), 
y  =  g(a,b,t)  to  get 


p  =  p'  (f(a,b,t)i  g(a,b,t)  ,  t) 

The  function  so  obtained  is  the  same  as  the  function  obtained 
frctn  the  Lagrange  approach. 

This  point  is  developed  in  some  detail  because  of  a 
certain  amount  of  confusion  in  the  literature  concerning 
Lagrangian  coordinates.  Many  statements  are  given  which, 
though  literally  accurate,  convey  a  misleading  impression. 

One  example  is  the  statement  that  the  Lagrangian  reference 
frame  is  not  inertial.  This  statement  is  literally  true, 
since  the  particles  do  not  change  their  initial  coordinates 
with  time,  regardless  of  the  pressures  acting.  In  that  sense 
the  coordinate  system  is  not  only  fixed  in  the  fluid  but  also 
distorting  as  well  as  moving,  with  the  fluid  and  hence  the 
Lagrangian  mesh  of  Figure  lb  is  not  inertial.  But  the  reader 
is  cautioned  against  the  error  of  thinking  that  the  frame  of 
reference  for  the  motion  is  moving.  The  frame  of  reference  is 
the  coordinate  system  with  respect  to  which  x  and  y  are  measured. 
These  are  the  Eulerian  coordinates  and  refer  to  a  coordinate 
system  fixed  in  space,  (the  fixed  axes  in  Figure  lb). 


While  the  scheme  of  using  the  initial  coordinates  for 
Lagrangian  variables  is  the  simplest  and  most  straightforward 
approach,  it  is  clear  that  any  function  of  the  initial  variables 
will  do  equally  well.  One  such  function  that  is  used  in  the 
scale  is  the  one  ve  will  describe  next. 
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Consider  the  equation 


M  = 


x(M,t) 

f  p(!>t)d| 

x(o,t) 


where  p  is  the  fluid  density  and  x,  the  coordinate  of  a  point 
in  the  fluid,  is  defined  as  a  function  of  M  hy  the  integral. 
If  we  take  a  partial  derivative  with  respect  to  M  on  both 
sides  we  obtan: 


1  =  p(x(M,t),t) 

We  consider  analogous  definitions  for  the  cases  of  a  sphere 
or  an  infinitely  long  cylinder  -  each  of  which  is  a  one  dimen¬ 
sional  problem. 


For  the  cylinder. 


R(M,t) 

M  =  2  rtf'  p(s,t)|d| 


and 


1=2*  p(R(M,t),t)  R(M,t)  |g 
For  the  sphere,  we  have 


r( M,t) 

M  =  4*  I'  p(|,t)|2d| 


1  =  4*  p(R(M,t),t  )  [R(M,t)]a  || 
where  R  is  the  appropriate  Eulerian  radius. 


We  now  introduce  a  quantity  designated  by  A  which  represents 
a  surface  area.  For  slab  symmetry,  we  choose  an  area  of  unit 
height  and  unit  width  (Ay  =  1,  Az  =  l)  and  obtain  A(M,t)  =  1. 

For  a  cylinder,  we  choose  unit  height  Az  and  angle  2jt  to 
obtain  a  surface  area:  A  =  2jt  R(M,t)  and  for  a  sphere  we 
choose  a  solid  angle  of  4jt  to  obtain: 


A  =  4jt[R(M,t)]2 


The  diagrams  are  as  follows: 


Cylinder : 


Sphere : 


A 


It  follows  from  the  law  of  conservation  of  mass  that  the 
mass  contained  inside  a  volume  bounded  by  the  same  fluid  particles 
is  constant  in  time  and  depends,  therefore,  only  on  the  original 
volume  and  the  coordinates  of  that  volume.  Hence  the  quantity 
M  may  be  chosen  as  a  Lagrangian  coordinate,  and  the  Eulerian 
coordinate  R  (or  x,  in  the  plane  case)  is  related  to  it  by  the 
differential  equation 
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3R  _  l 
Sm  p(M,t)  A(M,t) 


where  A(M,t)  =  1  for  a  slab,  A(M,b)  =  2-x  S(l4,t)  for  a  cylinder 
'and  4jt  R2  for  a  sphere. 

With  this  much  noted,  we  turn  our  attention  to  the  equations 
to  be  solved. 

II.  EQUATIONS  OF  THE  WUNDY  CODE* 

We  begin  our  discussion  of  the  Wundy  equations  with  some 
remarks  concerning  time  derivatives.  In  order  to  facilitate 
our  discussion  we  set  xx  =  x  x2  =  y  and  x3  ~  z  and  write  xt  as 
a  general  term  to  represent  Xj_,  x2  or  x3.  Similarly,  we  use 

a*  to  represent  ai,  a2,  or  a3  (ie;  a,  b  or  c). 

Now  consider  a  particle  located  at  the  point  (x.,,  x2,  x3). 

In  the  Lagrangian  formication  xj.  is  a  dependent  variable  and 
depends  on  aj  as  well  as  t.  But  the  at,  which  are  initial 
coordinates  (or  functions  thereof),  do  not  depend  on  time. 

Hence  we  write  x1l  =  and  xi  is  the  velocity  of  the  fluid 

St 

particle  at  the  time  t.  Similarly,  the  acceleration  is: 


Xi  = 


S^i 

St2 


In  the  Eulerian  formulation,  the  Xi  are  the  independent 
spatial  coordinates  but  depend  on  the  time.  The  velocity  and 
acceleration  are  now  given  by 

•  dx.  ,  ..  d2x. 

xjL  =  — —  and  Xi  =  i 

dt  dt2 


*The  reader  already  familiar  with 
skip  this  section. 


fluid  mechanics  may  easily 


the  difference  "between  the  two  forms  "being  the  same  as  the 
difference  "between  partial  and  total  derivatives  in  ordinary 
calculus. 

As  a  further  illustration,  consider  the  rate  at  which  the 
density  of  a  given  fluid  particle  changes  with  time.  In  the 
Lagrangian  formulation  this  is  simply  dp/dt.  In  the  Eulerian 
formulation  this  is1  1 ' _ 

n 

i 

dp  =  J.  V3 

dt  c^b  1=1  Sxj  dt 

or,  upon  replacing  dxi/dt  "by  the  Eulerian  velocities: 


&  =  k  +  £3  V, 

dt  St  i  ~  Sxj 

or  in  vector  form:  I  -%+V'o 

The  quantities  whose  time  derivatives  we  will  be  seeking 
are  usually  integrals  and  it  will  be  necessary  to  determine 
the  time  derivative  of  the  Jacobian  determinant  before  we  can 
handle  them.  We  suppose  the  Eulerian  coordinates  to  be  related 
to  the  Lagrangian  coordinates  by  a  set  of  equations: 
xA  =  Xj  (a^,  a2,  a3,  t)  and  the  Jacobian  determinant  of  these 
equations  is: 

J  =  Iti3I 

where  Tj. j  =  •  (^i/Saj)  and 

|  T±  j)  represents  the  determinant  formed  from  the  elements,  T^. 


1Note  that  the  x,y, z  co-ordinates  of  a  particle  are  themselves 
functions  of  time,  hence  the  total  derivative. 
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To  find  the  derivative  of  J,  we  write: 


dJ  _  3  dJ  a?rs 

^  "  \,s=1  ^7  dt 

The  partial  derivative  of  the  Jacobian  with  respect  to 
the  element  Trs  is  a  rather  complex  operation  and  will  be  set 
forth  in  Appendix  A  for  the  reader  familiar  with  tensor 
analysis. 

We  note  the  result: 


where  the  symbol  in  braces  denotes  the  rs  term  of  the  inverse 
of  the  matrix  formed  from  the  elements  Trs.  This  inverse  is 

simply  the  matrix  a .  Hence  for  our  time  derivative  we  have: 
3xr 


J23 

dt  r,s=l 


das  d  m 
dxr  dt  rs 


=  I3 

r,  s=l 


j  ^as  d 
c)xr  dt  5as 


=  r3 


3=1 


das  a  dxr 

oxr  das  dt 


- 


=  £ 


■S.. 


^as  5 
dxr  das  r 


gj 


With  this  derivative  out  of  the  way,  we  may  turn  our 
attention  to  the  evaluation  of  time  derivatives  of  integrals. 

To  this  end,  let  0  represent  any  physical  quantity  characteriz¬ 
ing  a  set  of  fluid  particles.  We  wish  to  integrate  0  over 
all  such  particles: 


1 = IS f  ^  dv 

c 

C,  being  the  region  comprised  of  such  particles.  As  the 
fluid  moves,  the  equation  of  the  boundary  changes  in  time  as 
does  0  itself.  But  the  Eulerian  coordinates  of  the  boundary 
are  functions  of  the  Lagrangian  coordinates,  a*,  and  we  can 
express  the  integral  as  an  integral  over  the  initial  volume  of 
fluid. 


1  =f If* 

c0 


JdV 


where  C0  is  the  boundary  at  time  zero  and  J  is  the  Jacobian 
of  the  transformation  from  the  at  to  the  x*  coordinates.* 

Quite  obviously, "J  depends  on  time,  but  the  bounding 
surface  no  longer  does  since  it  is  evaluated  at  t  =0. 


dl 

at 


0  J  dV 


*The  reader  unfamiliar  with  Jacobians  is  referred  to  Taylor: 
Advanced  Calculus,  pages  428-431. 
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Now  since  the  limits  of  integration  are  constants,  we  can  differ¬ 
entiate  under  the  integral  sign  (assuming  0  and  J  are  sufficiently 
smooth)  to  obtain 


(upon  inverting  the  transformation). 

This  provides  us  with  all  the  background  material  for  the 
conservation  laws. 

Conservation  of  mass: 
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and,  since  this  is  true  for  any  arbitrary  volume,  the  integrand 
must  be  zero. 


|L  +  (,  v.u  =  o 

Conservation  of  momentum: 


kill  p  dv  “  1  ^ 

c 

where  F^  is  the  k  component  of  the  force  acting,  and  the  sum 
is  over  all  forces. 


or  employing  the  mass  conservation  equations: 


c 


The  forces  acting  divide  into  two  kinds.  The  first  kind  is 
made  up  of  volume  forces  such  as  weight  and  electromagnetic  forces 
and  the  second  kind  is  made  up  of  forces  acting  on  the  volume  at 
the  surface.  The  hydrodynamic  approximation  consists  in  neglecting 
all  forces  except  those  arising  from  pressure. 

Since  the  pressures  are  on  the  order  of  megabars,  there  is 
no  problem  whatever  in  dropping  all  electromagnetic  forces  and 
weight  forces.  In  neglecting  all  viscous  forces  and  rigidity 
forces  in  the  metals,  the  program  is  perhaps  open  to  criticism. 
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but,  nonetheless,  the  approximation  is  not  as  bad  as  may  be  . 
anticipated  (Reference  l).  With  this  approximation,  we  write: 


III p  5T dv  "If  •pn,‘as 

C  Surface 

where  nk  is  the  k  component  of  the  normal  vector.  The  assumption 
of  the  hydrodynamic  approximation  has  the  effect  of  reducing  the 
stress  tensor  of  contimum  mechanics  to  a  diagonal  tensor  with  the 
elements  of  the  main  diagonal  all  equal  to  -  p,  (ie:  Ti j  =  -  P  8i j 
in  tensor  form). 

Applying  the  divergence  theorem  to  the  integral  at  hand, 
we  have: 


and  finally,  since  this  is  true  of  any  volume,  we  write: 


p  ^  +  Jr  =o 

St  Sxk 

which  is  the  form  of  the  momentum  equation  we  wish. 

The  energy  equation  is  also  easily  derived.  We  let  U  be 
the  internal  energy  per  unit  mass  and  let  h  be  a  vector 
representing  the  heat  current  out  of  the  surface.  (The  heat 
current  is  defined  as  that  vector  which  makes  the  integral 


h »n  ds 

Surf 
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equal  to  the  heat  that  is  transported  out  of  the  volume  through 
the  surface.)  Further,  suppose  that  heat  is  being  generated 
inside  the  surface  at  a  rate  X  per  unit  mass.  Then  applying 
the  first  Law  of  Thermodynamics,  the  rate  at  which  the  energy 
in  the  volume  changes  is  equal  to  the  work  done  on  the  body  by 
the  pressure  forces  plus  the  heat  generated  inside  the  body 
(X,  as  defined  above)  minus  the  heat  lost  through  the  surface 
(heat  current,  mentioned  earlier).  In  equation  form,  this  • 
statement  reads: 


_d 

dt 


J J J  p  (U  +  i  u2)  dv  =  -  J  J  p  (n  •  u)dA  -  J  J h  •  n  ds 


in 


+  111  p  X  dv 
c 

The  derivative  on  the  left  we  evaluate  as  before: 


III 


|p  (U  +  i  u2)  +  p  (U  +  uu)  +  7  •  u  p  (U  +  i  u2) 


dv 


-III 


p  (U  +  uu)  +  (U  +  i  u2)  (p  +  p  7 .  u)  dv 


-III 


p  (U  +  uu) 


dv 


where  the  mass  conservation  law  has  been  applied. 
Wow  we  have: 


J  J  J  p  (U  +  uu)dv  =  -  J  J  J  (Pu)  dv  “  III-  h  dv 

°  Vol  C  C 

'+///cX  dv 
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m 


p  (U  -  x)  +  puu  +  7  •  pu  +  7 •  h 


dv  =  0 


III 


p  (U  -  X)  +  u.  [pu  +  7p]  +  p  7  •  u  •+  7*  h 


dv  =  0 


+  p  7  •  u  +  7*  h 


dv  =  0 


where  the  momentum  conservation  law  has  been  employed.  Since 
the  integral  vanishes  for  every  volume,  it  follows  that: 


p  (U  -x)  +  p7*u  +  7*h  =  0 

For  the  Wundy  code,  the  heat  conducted  away  is  assumed 
negligible  and  h  =“  0.  The  source  of  heat  is  assumed  to  be 
viscosity  and  is  incorporated  into  the  pressure  term  by  the 
introduction  of  an  artificial  viscosity  term,  q  (Reference  2). 
This  term  is  the  very  core  of  the  solution.  By  manipulating  it 
properly,  we  arrive  at  a  continuous  solution  for  the  variables. 
The  shock  discontinuity  which  appears  is  smoothed  out  over  a 
few  zones  and  it  becomes  possible  to  find  difference  equations 
whose  solutions  will  approximate  the  revised  problem.  We-  will 
discuss  q  at  a  later  point  in  the  paper  but,  for  the  present,  we 
merely  note  that  -  p  X  =  q  7  •  u.  Solving  for  U,  we  obtain 
p  U  =  -  (p+q)  7  .  u.  Using  the  mass  conservation  equation,  we 
can  eliminate  7  •  u  to  get: 


U  =  (p+q)  -2  ^ 

P  St 

If  we  introduce  the  specific  volume  V'  =  1/p,  we  obtain: 


U  =  -  (p+q) 


SVl 

St 
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Summarizing,  our  equations  read: 


mass: 


momentum: 


energy: 


||+  p  7-u  -  0 

p*£-- JE 
St  Sx  k 


^  ~  “  <I*q) 


sv' 

sr 


The  first  of  these  equations  is  open  to  objection,  since 
there  is  little  point  in  obtaining  differential  equations  in 
the  density.  An  alternative  approach  is  to  introduce  the  mass 
variable  introduced  at  the  end  of  the  first  section.  Recall 
that 

x(M,t) 

M  =  J'  p(i*t)  A  ((•)  d| 

x(o,i) 


and 

Sx  _  1 
<3M  pA 

With  this  variable,  the  conservation  of  mass  is  automatic  since 
each  particle  is  permanently  identified  by  the  mass  contained 
between  the  initial  particle  and  the  particle  in  question. 

The  other  conservation  laws,  take  the  form: 

Momentum:  p^B=-^££!£pAor^ii=-A^E 
St  bx  SM  St  SM 
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To  obtain  the  equations  of  the  code,  we  must  multiply  this  last 
equation  by  pQ  and  set  pQU  =  E  and  p0V'  =  V.  Since  p0,  the 
initial  density,  is  constant  we  have: 


dE  _ 

St  " 


(p+q) 


dv 

St 


and  the  momentum  equation  is  unaltered: 


To  these  two  equations  of  motion  we  add  the  supplementary 
equations : 


V  =  p°/ p 
u  =  Su/St 

u  =  Sr/at 


Definitions :  Area 


Volume 


f  1 

slab 

< 

2  R 

cylinder 

1 

4  R2 

sphere 

r 

— 

R0  -  % 

slab 

«(H  f  -  R,*) 

cylinder  > 

|X3  -  Ki3) 

sphere 

(The  subscripts  i  and  o  refer  to  inner  and  outer  boundaries 
respectively  of  the  region  of  interest). 


Equation  of  State:  p  =  F(E,V) 


The  equation  of  state  is  a  thermodynamic  equation  relating 
(in  our  case)  the  pressure,  internal  energy  and  reduced  density. 
We  will  discuss  various  equations  of  state  somewhat  later  in 
this  report. 
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III.  DIFFERENCED  FORM  OF  THE  EQUATIONS 

The  basic  notion  behind  the  use  of  high  speed  computers  to 
solve  differential  equations  is  the  replacement  of  the  derivatives 
involved  by  difference  quotients  which  approximate  the  derivative. 
Thus  to  solve  the  equation: 


y(x0)  =  y0 


we  assign  A  x  a  value  and  compute  A  y  from  the  equation: 

Ay/Ax  =  x0.  The  variable  y  is  the  increased  from  y0  to  y0  +  A  y 
and  x  is  increased  to  x0  +  A  x.  The  process  is  then  repeated  and 
a  column  of  corresponding  x  and  y  values  are  produced. 

Similar  observations  apply  to  partial  differential  equations, 
except  that  there  are  more  variables.  To  this  end  we  number  the 
space  variable,  in  this  case  the  mass,  by  the  number  J.  The  time 

is  represented  by  the  number  n.  Thus  M-  is  represented  by: 

at 


30/dM  by:  ^  J+l  -  0  J  j  ^ (M. 


-  M,) 


To  facilitate  the  writing,  we  employ  0(J,N)  to  represent  0N 
hereafter. 


The  computation  proceeds  by  calculating  all  the  quantities 
desired  at  the  time  n  for  each  mass  cell.  The  spatial  computations 
are  performed  first  at  a  fixed  time,  then  the  time  is  advanced 
and  the  process  repeated. 

Fractional  indices  are  introduced  in  order  to  "center" 
the  difference  scheme.  In  this  fashion,  a  higher  accuracy  can 
be  achieved  with  the  same  number  of  zones.  The  quantity  0  (J  +  l/2) 
is  a  value  of  0  intermediate  between  0( J)  and  0( J+l) .  To  illustrate 
their  use,  let  us  take  a  differenced  form  of  the  equations  derived 
in  part  2.  We  write  these  equations  down  and  comment  on  them 
later. 
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Conservation  of  Energy: 


-  Cm) 


cSV 

5t 


•E(J  -  i/2,  fill)  =  E(J  -  l/2,  N) 


,lf 


T  i  P  (J  -  1/2,  N-)  +  P  (J 


Moment urn : 


4-  q(j  -  1/2,  N  +  1/2) 
-  V  (J  -  1/2,  N)j 


V  (J  -  1/2, 


N  +  1) 


JE  =  -  A 

cit 


DUDT  =  -  A(J,N)  Cpt<l).(  J+l/2,ri)  -  (p+q)(  J-l/2,N) 

m( J+l/2)  -  M( J-l/2) 


where  DUDT  is  the  computer  name  for  the  acceleration.  This 
expression  is  equivalent  to: 


DUDT  =  A(J,B) 

1/2  [m( j+1/2)  +  M( J-l/2)] 


Defining  Relations: 


V  =  p0/p 

V(J  -  1/2,  N  +  1)  =  po  Vol( J  -  1/2,  Nil)  /  m(J  -  1/2) 
_  c)x 

U  =  — 

dt 

x(J,  M  +  1)  =  x(j,  N)  +  u(J,  H  +  1/2)  .  A  Tmin 


IS 


-  1/2,  N  +  1) 


DUDT  =  ~ 

at 

u(J,  N  +  1)  =  u( J ,  N)  +  DUDT  A  Tnln 


1 

.Area  =  2jtX  A(J,  N  +  l)  = 

4jtX2 


Vol  = 


(Ro- 

*(  Ro2 


Ri) 

-  Ri2). 


4jt/3(R03  -  R^) 


1 

2jtX(«J>  N  +  1) 
4irX2(j,  N  +  1) 


,  X(J,  N  +  1)  -  X(J  -  1,  H  +  1) 

Vol  (J  -  1/2,  N  +  1)  =  *[  X2  (J,  N  +  1)  -  X2  (J  -  1,  N  +  1)] 

4jt/3  [X3  (J,  K  +  1)  -  X3  (J  -  1,  N  +  1)1 


These  are  the  basic  equations  of  Wundy  as  set  forth  in 
section  2.  The  superscript  n  refers  to  the  time  and  the  subscript 
j  to  the  distance.  The  subscript  j  -  i  identifies  a  variable  as 

in  the  space  bounded  by  the  two  lines  Rj  and  Rj_i . 

The  quantity  m(J  -  l/2)  is  the  mass  contained  in  the  zone 
between X(j)  andX(j  -  l).  It  is  obtained  from  the  previous 
mass  variable  M  by  the  relation  m(j  -  l/2)  *  M(j)  -  M(j  -  1). 

The  expression  M(j  +  l/2)  -  M(j  -  l/2)  should  be  taken  to  mean 
the  difference  between  two  masses.  The  first  mass  is  that 
contained  between  the  left  hand  boundary  of  the  initial  particle 
and  the  surface  corresponding  to  X(J  +  l/2).  The  second  is  the 
mass  between  the  same  initial  boundary  and  the  surface 
corresponding  to  x(j  -  l/2).  This  is  best  understood  by  reference 
to  the  following  diagram: 
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1/2) 


m(j  -  l/2)  m(j 


m 

M 

X(J  -  1)  X(J  -  1/2)  X< 

.J)  X(J  1/2)  X(J  +  1)  . 

[M(J  +  1/2)  -  M(J  -  1/2)] 


Since  M(  J  +  l/2)  -  M(j  -  l/2)  is  the  sum  of  the  masses  lying 
between  X(J  -  1/2)  and  X(j)  on  the  one  hand  and  between  X(j) 
and  x( T  -  1/2)  on  the  other.  If  x( J  -  l/2)  and  X( J  +  l/2)  ' 

divide  the  regions  X(j  -  1)  -  X(«J),  x(j)  -  X(j  +  1),  into 
two  halves  with  equal  masses,  then  the  mass  between  X(j  -  l/2) 
and  X(j)  will  be  half  the  mass  of  the  entire  region  x(J  -  1) 

-  X(«T),  that  is  l/2  m(j  -  l/2).  Similarly  the  other  mass 
involved  is  l/2  m(j  +  l/2),  so  that  M( J  +  l/2)  -  M(j  -  l/2) 

=  |  (m(J  +  1/2)  +  m(j  -  1/2)). 

With  this  in  mind,  it  is  obvious-  that  the  two  expressions 
for  DUDT  in  the  difference  equations  are  equivalent. 

At  this  point,  the  question  arises  whether  the  differential 
equations  could  be  represented  by  another  set  of  difference 
equations.  The  answer  to  this  question  is  yes.  Fromm  (Refer¬ 
ence  8)  considers  several  schemes  of  differencing  and  evaluates 
them.  Unfortunately,  Fromm's  criterion  for  evaluation  is 
basically  a  stability  criterion.  He  considers  the  best  equation 
to  be  the  one  that  allows  the  largest  A  t  to  be  used  before 
instability  sets  in.  While  this  may  be  a  good  mathematical 
criterion,  the  obvious  physical  criterion  of  agreement  between 
calculated  and  observed  values  has  not  yet  been  fully  employed. 
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As  an  example  of  the  alternatives  available,  the  following 
extract  from  Fromm's  work  is  given: 

1)  X(J  -  1/2,  N  +  1)  =  X(J  -  1/2,  N)  +  u(j  -  1/2,  H)  At 

2)  X(J  -  1/2,  N  +  1)  =  X(J  -  1/2,  N) 

+  |  [u(j  -  1/2,  N  +  1)  +  u(j  -  l/2,  N)]  At 

3)  X(J  -  1/2,  N  +  l)  =  X(j  -  1/2,  N)  +  u(j  -  1/2,  N  +  1)  At 

Of  these  equations,  Fromm's  work  indicates  that  the  last  is 
preferable.  On  physical  grounds,  however,  it  seems  that  some 
sort  of  average  velocity  would  yield  the  best  results  in  the 
computation  of  the  final  position,  since  this  would  represent 
the  velocity  of  the  zone  at-  some  intermediate  time.  There  is 
also  a  tendency  to  choose  the  initial  value  of  velocity  partly 
because  we  ordinarily  extrapolate  from  present  knowledge  (u(n) 
and  X(N))  to  future  knowledge  X(N  +  1).  Physically  Fromm's 
choice  of  the  best  equation  runs  counter  to  the  normal  choice. 


By  considering  a  Taylor  expansion  of  X  as  a  function  of 
time,  we  will  gain  seme  insight  into  the  difference  equations. 
Accordingly,  we  write:1  (using  superscripts  and  subscripts 
throughout  for  clarity. 


*+.l/  2  .  -v  N+V2  P  o 

.At  _t .  At2  dsx 

,.i/s  2  F/,-1/1  ~  |st*j 


H  N+1/2  IN+1/2  „  „  . N+1/ 2 

X  =  X  _  A£  dx  +  t2  d2X  1  + 

,.1/s  ,.1/a  2  St|,-l/a  B-  a^),.l/a  +  ••• 

Upon  subtracting  these,  we  obtain: 

N+1/2  N  N+1/2 

x  .  "  x  =  At  ~  ,  +  0(At)3  and  the  most 

j-1/2  j-1/2  &  j-l/2 

accurate  representation  for  fixed  At  ought  to  be: 

xNote  that  fundamental  time  increment  At  is  the  time -from  cycle 
N  to  cycle  N+l.  Since  we  have  used  cycle  N+l/2,  the  time 
difference  between  the  2  cycles  is  At/2. 
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H+l 


II 

X  +  At 

J-l/2 


H+1/2 

J-1/2 


N  N+-/2 

=  X  +  At»u  . 

J-l/2  J.I/2 


But,  as  Fromm  points  out,  the  error  introduced  by  neglecting 
terms  of  higher  order  may  actually  be  lower  in  the  case  where 
uI,+1  because  the  neglected  terms  tend  to  cancel  each  other  out 
in  one  case  and  not  in  the  other. 

Similar  differencing  schemes  may  be  tried  on  all  the 
equations,  but  the  scheme  used  in  Wundy  is  the  only  one  that 
need  concern  us  here. 


Since  the  difference  equations  are  not  the  same  as  the 
differential  equations  they  represent,  but  only  an  approximation 
to  the  differential  equations,  the  solutions  to  the  difference 
equations  are  not  the  same  as  the  true  solutions  but  only  an 
approximation.  It  is  reasonable  to  suppose  that  the  solutions 
to  the  difference  equations  can  be  made  as  close  as  we  please 
to  the  true  solutions  by  choosing  At  sufficiently  small  and 
making  m(J  +  l/2)  sufficiently  small.  But  non-linear  differ¬ 
ential  equations  are  rather  ornery  oujeets  and  have  a  way  of 
confounding  predictions  of  this  kind.  In  our  case,  however,  the 
two  solutions  can  be  made  arbitrarily  close  providing  only  that 
(l)  At  and  m(j  +  l/2)  are  chosen  sufficiently  small  and'  (2) 
that  the  system  of  difference  equations  is  stable,  the  theorem 
(stating  that  stability  of  a  system  of  difference  equations  is 
necessary  and  sufficient  for  convergence  of  the  solution  of 
these  equations  to  the  true  solution)  being  due  to  Peter  Lax 
(References  3  and  6). 

The  stability  of  a  system  of  difference  equations  refers 
to  the  amplification  of  errors  introduced  into  the  computing 
process  by  various  means  -  round-off  errors,  inaccurate  data 
and  so  on.  If  such  errors  are  amplified  and  grow  in  time  as  the 
computation  proceeds,  the  difference  equations  are  said  to  be 
unstable.  If  not  the  equations  are  stable.  Striking  examples 
of  numerical  instability  are  given  by  Richtmeyer  (Reference  3) 
and  the  following  example  is  taken  from  the  source. 
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The  continuous  lines  represent  the  solution  to  a  differ¬ 
ential  equation  and  the  dashed  lines  the  solution  to  the 
difference  equations.  It  is  obvious  that  the  two  solutions 
are  not  approaching  one  another  as  t  increases. 

Instability  is  avoided  by  insuring  that  a  certain  relation 
holds  between  the  increments  to  the  distance  and  time  variables. 

C our ant,  Friedrichs  and  Lewy  deduced  the  first  stability  cri¬ 
terion  and  Von  Neumann  and  Richtmeyer  have  applied  this  criterion 
to  the  Wundy  code.  The  criterion  is  automatically  satisfied 
by  the  DELTA!  subroutine  which  computes  At  -  the  change  in 
time  from  1  cycle  to  the  next.  It  is  important  to  note  that 
no  stable  set  of  difference  equations  exists  to  describe  the 
problem  with  a  shock  discontinuity  present.  This  is  why  the 
artificial  viscosity  was  introduced. 

The  functional  form  of  the  artificial  viscosity  is  rather 
arbitrary  but  the  following  conditions  are  imposed: 

(1)  When  the  artificial  viscosity  term  is  incorporated 
into  the  problem,  the  equations  must  have  continuous  solutions. 

(2)  The  distance  over  which  the  shock  is  spread  must  be 
everywhere  on  the  order  of  the  thickness  of  the  zones  used  in 
the  numerical  computation,  independently  of  the  shock  strength 

( so  that  the  same  q  will  work  for  all  shocks)  and  independently 
of  the  condition  of  the  material  into  which  the  shock  is  moving 
(so  that  q  is  determined  solely  by  the  variables  which  influence 
the  pressure,  density,  velocity  and  the  other  physical  quantities 
Involved) . 

(3)  Outside  the  shock  front,  the  effect  of  q  must  be 
negligible.  This  is  to  be  taken  to  mean  that  the  difference 
between  the  quantities  computed  with  the  use  of  q  and  the 
quantities  computed  by  more  exotic  schemes  will  be  negligible, 
which  implies  that  q  itself  is  small. 
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(4)  The  Hugoniot  equations  must  hold  when  all  other 
dimensions  are  large  compared  to  the  shock  thickness. 

von  Neumann  and  Richtmeyer  discovered  that,  for  1 
dimensional  flow  in  a  substance  which  obeyed  the  y  law 
expression: 


<1 


(CAx)2  du 
y '  Sx 


satisfied  all  these  conditions.  A  later  correction  was  made 

setting  q  =  0  for  ^  s  0.  For  if  <Wdx  is  non-negative  then 
dx 

the  front  of  the  zone  is  moving  more  rapidly  than  the  rear 
and  the  zone  is  expanding.  In  this  case,  no  shock  is  formed 
or  likely  to  form  and  consequently  the  undisturbed  equations 
have  a  continuous  solution  and  there  is  no  need  of  the  smoothing 
effect  of  q.  The  form  of  q  used  in  Wundy  is: 


q  =  -  §u 

.  H  !5x'  Sx 


But  V'  =  v/p0  so  that 


if  <  0 

dx 


4-  fcsEp. 


du 
l  dx 


q(j  -  1/2,  N  +  l/o) 

=  a  m  Cg  [u(j,  N  +  1/2)  -  u(J  -  1,  N  +  l/g)]2 
0  V(J  -  1/2,  N  +  1/2) 


=  C2  P°  (I)  [u(J,  N  +  1/2)  -  u(j  -  1,  N  +  l/2) ] 2 
|  [V(J  -  1/2,  N  +  1)  +  V(J  -  1/2,  N) 


Using  the  code  symbol  CQSQX4,  we  set  von  Neumann's  C2  =  CQSQX4/4 
to  obtain: 
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q(j  -  l/2,  N  +  1/2) 


_  CQSQX4  Pp  (l)[u(j,  N  +  1/2)  -  u(j  -  1,  H  +  l/2)]; 
4  |rv(J-l/2,  N  +  l)  +  V(J  -  1/2,  N)] 


=  CQSQX4  Po  (l)[u(j,  N  +  1/2)  -  uf  J  -  1,  N  +  l/2)]2 
~  2  “  WTj  -1/2,  N  +  1)  +'  VI J  -  1/2,  N)] 

and  CQSQX4  (i)  is  set  equal  to  9  for  explosives  and  16  for  metals 
in  the  present  code. 

IV.  EQUATIONS  OF  STATE 

The  equations  of  state  are  normally  determined  empirically. 
For  explosive  gases,  the  equation  used  is  the  ideal  gas  equation: 


pV  =  RT;  E  -  E0  =  CVT 

where  E  is  the  internal  energy  per  mole  and  p,V,T  R  and  Cv  all 
have  their  usual  significance  and  refer  to  1  mole  of  gas. 
(Refer  to  Zemansky:  Hera,  and  Thermodynamics,  page  120)  Recall 
Cp  -  Cv  =  R  and  CP/.C.  =  .  and  we  have: 


Dividing  "both  .'  ides  hy  the  atonic  weight  in  grams,  on  the  left 
we  obtain  the  al  energy  per  gram  instead  of  per  mole) 

and  on  the  right,  we  hsve  the  volume  of  one  mole  divided  by 
the  mass  of  one  mole  wnich  is  the  reciprocal  of  the  density. 

0  -  u0  =  ;-rr  p/p 
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Multiply  by  p0  -  the  initial  density  and  the  result  is: 

E  -  Ec  =  p0  (U  -  U0)  =-1-5  (p  /p) 

7-1 


E  -  E0  =  — ■ i—  pV 
7-1 

where  V  =  Po/p  as  defined  earlier. 

This  equation  represents  the  energy  of  an  ideal  gas  for 
7  =1.4.  For  values,  of  y  in  the  neighborhood  of  3*,  the  equation 
also  represents  not  only  a  condensed  explosive  but  also,  to  a 
good  approximation  the  high  pressure  gases  released  by  the 
detonation  process  (Reference  7).  The  energy  that  is  released 
into  the  wave  is  controlled  by  the  burn  fraction,  F,  which  is 
defined  as  the  following  ratio: 


where  V  =  Po/p  and  Vcj  =  P°/pcj  the  subscript  CJ  referring 
to  the  Chapman- Jouget  conditions  (Reference  4).  If  the 
Chapman- Jouget  hypothesis  is  satisfied,  then: 

(po/p)  =  — L_  =  vCJ 

7  +  1 


and 


1 

7  +  1 


p  =  Az  l  =  (7  +  1)  (1  -  V) 

1  -  VCJ 

which  is  in  the  form  used  in  the  code.  The  fraction  of  the  energy 
released  is  then  F. 


*This  value  of  7  need  no  longer  represent  the  ratio  of  specific 
heats  in  the  explosive 


In  the  actual  working  of  the  code,  the  "burn  fraction  is 
subjected  to  several  tests  designed  to  insure  that  the  computed 
situation  is  physically  realistic.  If  F  <  .00001,  it  is  assumed 
that  the  difference  between  F  and  0  is  due  to  errors  in  the 
calculation.  If  F  >  1,  it  is  assumed  that  the  difference  in 
F  and  1  is  due  to  the  same  cause.  F  is  then  set  equal  to  1. 

If  F  begins  to  decrease  in  time,  F  is  set  equal  to  1.  Finally 
if  F  is  equal  to  1  in  every  zone,  the  energy  of  the  explosive 
is  spent  and  the  routine  is  skipped  on  the  next  cycle. 

The  equation  of  state  and  the  energy  equation  are  both 
used  in  advancing  the  energy.  From  the  above  considerations,  we 
have: 

3>p/dE  =  (7  -  1)/V 
or  in  differenced  form: 

DPDE  =  (7  -  1)/V  (J  -  1/2,  N  +  1) 

The  energy  equation  of  section  II  is  now  employed 
dE/cit  =  -  (p  +  q)  dV/dt 

to  obtain: 

E(J  -  1/2,  N  +  1/2)  -  E(J  -  1/2,  N) 

=  -  (p  +  q)(J  -  1/2,  N)  [V(J  -  1/2,  N  +  1) 

-  V(J  -  1/2,  N)] 

But,  inside  the  explosive,  this  equation  is  subject  to  error 
because  energy  source  terms  are  not  included.  Moreover  the 
pressure  p(j  -  l/2,  N)  is  decidedly  not  the  pressure  acting 
over  the  entire  time  cycle  or  even  a  good  approximation  to 
it  since  the  pressure  can  be  expected  to  increase  greatly 
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when  the  energy  of  the  detonation  is  released.  For  these  reasons, 
an  iteration  process  is  used  in  computing  E(j  -  l/2,  N  +  l). 

First  we  set 


El  =  E(J  -  1/2,  N)  -  (p  +  q)(j  -  1/2,  N  +  1)  [V(J  -  l/2,  N  +  1) 

-  V(J  -  1/2,  M) ] 

PI  =  DPDE  •  El  •  F(J  -  1/2,  N  +  1) 

These  intermediate  pressures  and  energies  are  then  employed 
in  obtaining  the  resultant  pressure  and  energies  as  follows: 

E(J  -  1/2,  N  +  1)  =  E(J  -  1/2,  N)  -  |  [PI  -  P(J  -  1/2,  N)  ] 

•  [V(J  -  1/2,  N  +  1)  -  V(J  -  1/2,  N)] 

P(J  -  1/2,  N  +  1)  =  DPDE  •  E(J  -  1/2,  N  +  1)  *  F(J  -  l/2,  N  +  1) 

The  computation,  having  advanced  E  and  P,  returns  to  the 
main  program. 

Of  all  the  inert  materials  a  genuine  equation  of  state 
is  available  only  for  aluminum  as  of  this  writing.  For  other 
materials,  the  Hugoniot  equation  is  used.  This  equation  is 
strictly  an  equation  of  process  and  not  of  state,  the  process 
involved  being  a  shock  transition.  The  Hugoniot  is  accurate, 
therefore,  only  for  shock  transitions,  but  very  weak  shocks  are 
nearly  isentropic  and,  if  the  drop  in  pressure  from  zone  to  zone 
is  not  too  great  in  the  relief  wave,  the  Hugoniot  will  be  a  good 
approximation  to  the  adiabat  which  is  the  actual  curve  along  which 
the  relief  transitions  take  place. 

The  Hugoniot  (Reference  5)  may  be  obtained  from  a 
straightforward  empirical  c urve  fitting  process  or  by  utilizing 
the  empirical  observation  that  the  shock  velocity  is  a  linear 
function  of  the  particle  velocity  at  sufficiently  high  pressures. 
The  Rankine-Hugoniot  equations  can  then  be  used  to  determine 
the  equation  of  the  Hugoniot  curve  in  P  -  q  coordinates  (where 
H  is  the  quantity  (p/p0  -  1)).  A  typical  calculation  of  this 
kind  is  presented  in  Appendix  B. 
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The  pressure  and  internal  energy  are  now  advanced  in  a 
manner  similar  to  that  employed  for  explosives.  The  intermediate 
Variables,  p(p)  (in  the  code  language  POFMU)  and  p  (or  DMU)  are 
introduced  and  the  empirical  relation  between  them: 


POFMU  =  (C  +  C2  DMU)  DMU 

El  =  E(J  -  l/2,  N)  -  [P(J  -  1/2,  N)  +  q(j  -  1/2,  N)] 

*  [V(J  -  1/2,  N  +  1)  -  V(J  -  1/2,  N) ] 


E(J  -  1/2,  N  +  1)  =  El  -  i  [POFMU  +  CpEI  -  P(j  -  l/2,  M)  ] 

2 

.  [V(J  -  1/2,  M  +  1)  -  V(J  -  1/2,  N)] 

P(J  -  1/2,  N  +  1)  =  POFMU  +  C2  E(J  -  1/2,  N  +  1) 

If  the  computed  pressure  turns  out  to  be  less  than  -3  kilobars, 
it  is  set  equal  to  -3  kilobars.  Since  a  tension  (negative 
pressure)  of  3  kilobars  vail  cause  many  metals  to  fracture, 
a  greater  tension  cannot  appear  in  the  physical  explosion. 

This  introduces  the  negative  pressure  cutoff. 

V.  THE  PROGRAM 


The  computation  begins  by  reading  off  the  initial  values 
of  the  variables  from  a  set  of  IBM  cards.  The  input  variables 
are: 

(1)  Equation  of  state  for  region  J. 

(2)  Initial  density  of  region  J. 

(3)  The  constant  CQSQX4  for  region  J. 

(4)  The  constand  7  for  region  J. 


30 


(5)  Certain  control  variables,  notably: 

(1)  MURIN 

(2)  mv 

(3)  NALTQ 

(4)  NSWEEP 

(5)  K 

This  list  of  control  variables  is  by  no  means  exhaustive. 
These  are  given  as  a  sample.  The  control  MURIN  defines  the 
initial  conditions  at  the  innermost  boundary.  If  the  geometry 
is  cylindrical  or  spherical,  MURIN  =  1.  For  rectangular  (slab) 
geometry,  there  are  two  possibilities,  either  a  rigid  wall  at 
the  origin  (set  MURIN  =  1)  or  a  free  surface  MURIN  =  0.  INV 
is  a  control  that  determines  the  initial  state  of  the  material. 
If  the  material  is  in  its  normal  state,  INV  =  0.  If  it  is 
desired  to  have  the  material  in  a  compressed  state  when  the 
program  starts,  INV  is  set  equal  to  any  number  other  than 
zero,  (for  example  IN'/  =  1).  Then  the  ratio  of  the  density 
in  the  compressed  state  to  the  uncompressed  density  is  read  in. 

NALTQ  is  set  equal  to  zero  if  the  von  Neumann  form  of 
the  artificial  viscosity  q,  is  desired.  If  NALTQ  ^  0,  then 
a  linear  q  will  be  used.  This  form  of  q  has  been  used  in 
certain  underwater  explosion  problems  at  Naval  Ordnance 
Laboratory,  White  Oak.  NSWEEP  is  a  control  which  introduces 
a  limited  calculation  procedure.  NSWEEP  =  0  causes  the 
velocities  of  three  adjacent  zones  to  be  tested.  If  the  sum 
of  the  absolute  values  of  these  three  zones  is  less  than  1CT6, 
and  the  time  cycle  is  less  than  2,  then  it  is  assumed  that 
the  shock  has  not  yet  reached  those  zones  and  the  computation 
is  stopped  at  the  last  of  these  zones. 

The  last  control  mentioned,  K,  determines  the  geometry 
of  the  problem.  K  =  1  for  rectangular,  2  for  cylindrical, 
and  3  for  spherical  problems. 

When  all  the  input  values  have  been  read  in,  all  the 
variables  of  the  calculation  (N,  J,  J  +  l/2,  etc.)  are  set 
equal  to  their  initial  values  and  the  patch  subroutine 
(PCH),  is  called.  This  subroutine  is  employed  to  begin  the 
process.  The  command  PCH(lJ)  =  TYME  alters  the  value  of  the 
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variable  in  the  location  IJ  to  a  new  value  designated  by  the  code 
word  TYME.  It  has  been  found  most  convenient  to  change  the  velocity 
of  the  first  interface  to  a  suitably  high  value  in  order  to  begin 
the  detonation  process.  With  the  first  zone  being  compressed,  the 
burn  fraction  rises  and  energy  is  released.  This  builds  up  the 
pressure  and  the  detonation  process  is  underway.  Next,  the 
subroutine,  GENSUB,  is  called. 

This  subroutine  sets  a  series  of  control  variables  which 
locate  the  various  interfaces  and  any  voids  that  may  be  present. 

The  program  is  now  ready  to  go.  The  limited  computation  test 
is  applied  to  determine  an  upper  limit  for  the  calculations. 

MURIN  is  then  tested  to  determine  the  fate  of  Hie  first  zone. 

For  a  rigid  wall  condition,  u(J,  N  +  l/2)  =  u( J,  N  -  .l/2) 

=  TYME.  For  a  free  surface,  EUBT  is  computed  and  u  and  x  are 
advanced.  The  computation  now  proceeds  to  advance  each  zone 
in  accordance  with  the  difference  equations  cited  in  section 
III.  The  velocity,  coordinate,  area,  reduced  density  and 
artificial  viscosity  are  all  computed  for  each  zone.  At  the 
void  closures,  the  special  routines  discussed  in  section  VI 
are  employed.  When  these  computations  have  all  been  performed, 
the  equation  of  state  is  selected  and  the  pressure  and  internal 
energy  are  both  advanced.  The  subroutine  DELTAT  is  called  up 
and  a  new  value  of  At  is  generated.  This  value  is  calculated 
from  stability  considerations. 

The  time  is  now  advanced  from  t  to  t  +  At  and  the  energy 
is  checked  to  insure  that  the  total  energy  does  not  depart 
from  the  initial  value  of  the  energy  by  more  than  10$.  If  it 
does,  the  computation  is  stopped  and  E  WRONG  is  printed  out. 

(This  energy  check  is  optional  and  is  controlled  by  NCHEKE 
which  is  set  equal  to  1  if  this  check  is  desired,  and  0  if 
the  check  is  to  be  bypassed.)  The  computer  now  compares 
t  +  At  with  the  time  limit  set  by  the  operator.  If  the  run 
time  is  larger  than  t  +  At,  the  computation  is  started  all 
over  again  using  this  advanced  value  of  t.  If  not,  the 
output  is  summoned  and  the  results  of  the  computation  are 
printed  out, 

VI.  THE  VOID  SUBROUTINE 

The  calculations  outlined  thus  far  proceed  very  well  in 
normal  situations,  but  occasionally  an  extraordinary  situation 
arises  and  special  attention  is  required.  The  commonest  of 
these  situations  occurs  when  a  void  is  enclosed  between  two 
regions. 
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A  new  Interface  must  be  introduced  (XVOID)  and  appropriate 
corrections  must  be  made  in  the  equations.  The  interface  on  the 
right  ifcere  the  material  begins  again  is  subscripted  JVOID  and 
the  two  collide  after  a  time.  The  geometry  is  as  follows: 


1 

Region 

VOID 

Region  I  +  1 

I 

Origin  XVOID  X( JVOID, N) 


The  interface  XVOID  is  advanced  as  follows: 

VD  DUDT  =  :  (p4q)(JV-l) 

VD  AREA 

|  m(  J  -  1/2) 

UVCID  (I,  N  +  1/2)  =  UVOID  (i,  N  -  l/2)  +  AT  VDDUDT 

XVOID  (I,  N  +  1)  =  XVOID  (I,  H)  +  UVOID  (l,  N  +  l/2) 

The  first  of  these  equations  corresponds  to  the  ordinary 
differenced  form  of: 


du  =  - 

at 

except  that,  for  the  void,  the  pressure  in  zone  J  -  l/2  is 
labelled  p(JV  -  l)  and  the  pressure  in  zone  J  +  l/2  is  set 
equal  to  zero  (since  there  is  no  pressure  in  a  void). 

Similarly  the  expression  for  mass  is  altered  since  the  zone 
mass  of  the  void  is  zero.  Thus,  only  m(j  -  1/2)  appears. 

The  right  hand  boundary  is  advanced  in  a  similar  manner: 
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u(JV,  N  4.  1/2)  =  u(JV,  M-1/2)  ~  VJDUDT 


X(jV,  N  +  1)  =  X(JV,  N)+£T  u(JV,  N  +  l/2) 

The  two  X's  are  now  compared  to  investigate  whether  the  material 
on  the  left  has  overtaken  the  material  on  the  right  or  not.  If 
XVOID  is  still  smaller  than  X( JVOID),  then  the  void  is  still 
open  and  the  computation  proceeds  hy  advancing  X(JV  -  1,  N  +  l), 
computing  a  new  value  of  q  and  returning  to  the  main  program. 

(X(JV  -  1,  i!  fl)  is  advanced  here  in  order  to  compute  q. 

It  is  recomputed  later,  hut  it  was  though  best  to  keep  all 
the  void  calculations  together.) 

If  XVOID  >  X(JV,  N  +  l),  a  void  closure  has  occurred  and 
the  program  proceeds  to  compute  a  new  time  increment  and 
employs  this  in  subsequent  calculations.  XVOID  is  now  set 
equal  to  zero  and  the  void  index  JVOID  is  set  equal  to  -JVOID 
which  serves  as  a  test  for  void  closure.  Hereinafter,  the 
subroutine  is  skipped. 

UVOID  is  advanced  from  the  VDDUDT  computation  and  U(JV,  N  +  l/2) 
is  also  advanced.  After  X(<JV,  N  +  l/2),  the  common  interface,  is 
advanced,  a  quantity  called  UTEMP  is  introduced.  This  is  the 
common  velocity  of  the  two  zones  adjacent  to  the  interface 
X(JV,  N  +  1)  as  computed  from  the  usual  formula  for  inelastic 
collision: 


1  1 
UTEMP  =  ?  m(jV)  U(JV,  N  +  l/2)  +  ?  m(«JV  -  l)  UVOID 

|  m(  JV)  +  |  m(JV  -  1) 

X(JV  -  1,  N  +  1)  is  now  advanced  just  as  before  and  areas, 
volumes  and  reduced  densities  are  canputed  and  the  computation 
is  returned  to  the  main  program. 

In  the  main  program,  itself,  a  further  correction  is  made 
to  adjust  the  energies  of  these  two  zones.  The  mechanical 
(kinetic)  energy  lost  by  the  first  zone  is: 


|  m(JV  -  1)  (VOID2  -  UTEMP2) 
and  that  gained  by  the  second  zone  is: 


|  m(  JV)  (IJTEMP^-  u  ( JV,  N  +  l/2)2) 
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The  mechanical  energy  lost  is  the  difference 


|jm(JV  -  1)  ( UVOID2  -  UTEMP2)  +  m(JV)  (U(JV,  N  +  l/2)2  -  UTEMP2)! 

This  will  appear  as  thermodynamic  energy  and  a  plausible  assumption 
for  metals  is  to  assume  that  this  energy  is  divided  equally 
between  the  two  zones  concerned.  If  we  neglect  the  slight  change 
in  density  of  the  metals,  we  have: 

m(  JV)  (E'  -  E2yp2  =  Total  change  in  internal  energy  of 

zone  JV.  x  mass  /  density) 


K  "  E2=  (paM^V)  l/e{m(JV  -  1)  [UVOID2  -  UTEMP2] 
+  ra(jV)  [U(JV,  N  +  l/2)2  -  UTEMP2]! 


‘4  {  I  P2  N  +  Vs)2  -  UTEMP2]  +  |  p2  [ UVOID2-UTEMP2 ]^- 


And  similarly 


Ei  -  El  = 


“ —  —Tm( 
-1)  8^ 


m(jV-l)  8 


JV  -  1)  [UVOID2  -  UTEMP2] 

+  m(<JV)[U(JV,  K  +  1/2)2  -  UTEMP2]! 


=i  ||  P1  (UVOID2  -  UTEMP2)  +  |  Pi  (u(  JV,  N  +  l/2)2  -  UTEMP2)! 


For  a  gas  closing  on  a  metal,  the  low  thermal  conductivity  of 
the  gas  is  presumed  to  prevent  any  appreciable  transfer  of  energy 
to  the  metal  so  that  we  may  neglect  £E2  and  employ  the  second 
equation  given  above  together  with  the  approximation  UVOID  »  UTEMP 
and  UVOID  »  U(JV,  N  +  l/2).  Neglecting  both  these  factors  and 
doubling  the  energy  on  the  right  in  order  to  insure  that  all  the 
mechanical  energy  lost  is  converted  into  thermal  energy  of  the 
gas,  we  obtain: 
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E,'  =  E,  +  -  -  -  P-a  -  [m(  JV  -  1)  UVOID2] 
4  m(«JV-l) 

=  Ea  +  -5  (1/2  pg)  UVOID2 


•which  is  the  equation  used  in  Wundy. 

When  these  adjustments  have  "been  made,  JVOID  is  set  equal 
to  zero  and  the  adjustment  is  skipped  thereafter. 

VII.  QUESTIONS  ON  QUASI-WUNDY 


There  are  several  questions  that  arise  in  the  course  of 
working  on  the  Quasi-Wundy  code  and  these  will  he  considered 
here. 


The  first  of  these  questions  concerns  the  necessity  for 
mass  matching  at  interfaces  between  adjacent  zones.  This 
question  ordinarily  arises  when  a  very  light  metal  is  placed 
in  contact  with  a  very  heavy  one  and  the  number  of  zones 
required  in  order  to  obtain  a  good  notion  of  the  processes 
occurring  in  the  light  metal  makes  the  number  of  zones  in  the 
heavy  metal  prohibitively  expensive  if  the  masses  of  adjacent 
zones  are  matched.  As  an  example,  consider  a  thin  aluminum 
slice  -  2.5xl0“3  cm  thick  -  in  contact  with  a  copper  plate. 

The  density  of  aluminum  is  2.785  gm/cm3  and  that  of  copper 
is  8.93  gm/cm3.  It  is  estimated  that  not  less  than  10  zones 
are  needed  in  the  aluminum  in  order  to  get  a  reliable  view  of 
the  motions  of  the  shock  waves  through  the  aluminum.  The 
thickness  of  each  zone  of  aluminum  is  2 .5x10" 4  cm.  Let  us 
consider  a  unit  cross  section  and  the  mass  of:  each  zone  is 
2.5xl0~4  x  1  x  2.785  =  6.9625  x  10”4  gm.  The  copper  plate 
is  .635  cm  thick  and  if  we  require  that  the  mass  per  zone  be 
6.9625  x  10" 4  gms,  we  have,  for  the  number  of  zones: 

N  =  Total-  mass  of  Copper  =  1  x  .655  x  8.95  x 
mass  of  zone  6.9625 

N=“  8,000 

The  cost  of  computation  for  8,000  zones  is  rather  high.  Thus  the 
question  naturally  arises:  Is  it  necessary  to  insure  that  the 
masses  of  the  zones  are  equal  or  nearly  so? 
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Consider  the  Taylor  expansion  of  the  pressure  near  a  point, 
m0,  at  which  the  increment  Am  is  changed.  (Normally  m0  will  he 
the  interface  between  two  regions).  According  to  a  well  known 
formula: 


p(m)  =  p(m0)  +  (&)  +i_(^)2(^!|)  +i_(M)3(^£)  +  ... 

*  °'  2  W»o  21  V2  '  V&n2'»o  31  3  '  ^m3  0 


where  m  =  m0  ±  l/2  Am 

Designate  (m0  +  i  Ag  m)  by  P  J+l^  p(m0  -  i  &:  1)  by  p  O.  and 

d.  2  £  —  "~2 

(m0)  =Pj.  Then  we  have: 


p,,=pj  +  fa  2£+i_  (<%t  (S£,a  +  1.  (£&  (2£,a+.. 

J  WJ  2  2:  dnrJ  V  2  '  3  ^m3'*  2  ' 


P  x  =  P  - 
0>H-i 


££  +  i_  (&  {t^ .  i.  (&  (SA)3  +  ... 

dm J  2  2i  dm2  J  2  31  5m3  J  2 


Subtracting  the  first  expression  from  the  second: 


P  j-irp  j4. 

2  2 


'  Axm  +  Agm  1 
2 


-  —  (III)  {  (Aim)3  +  (A?m)3  1  + 

48  ^clm3yJ  l  1  2  i 


Divide  both  sides  by  i  (A^m  +  Agin) 


and  we  have: 


Pj4 


+P, 


4 


~  (Axm  +  Agin) 


=  (fe  1 

dm'j  4 


(Aim  -  A2m) 
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[(Aam)2  -  (Aam)  (Agin)  +  (Agm)2]^-... 


If  A^m  =-  A .  tn,  A^m  -  A 2m  ■=-  0  and  our  equation  reads: 


1) 


pj4  - 


—  -  <J&*  *  ST  <&»  0t“2) 


|  (Axm  +  Agm) 


!SI 

dm5 


The  approximation  used  in  the  code  is  then  simply  to  drop  terms 
of  order  greater  than  the  second  and  set 


2) 


pj-J  “^4 


-  ($E)  _ 

dm  j  1 

-  (Axm+  A2m) 


and  employ  this  in  the  equation 


ou  ^  _  6  di 
<§t 


=  -  AfE 

.  dm 


which  advances  the  velocity  and  then  the  X  coordinate  and  hence 
the  density  and  indirectly  the  pressure  in  the  next  zone. 

If  the  masses  are  not  matched,  this  equation  can  he 
seriously  in  error. 

In  an  actual  run  (Quasi-Wundy  84),  an  attempt  was  made  to 
match  the  masses  of  the  zones  in  the  aluminum  foil  mentioned 
above  and  the  copper  plate.  In  order  to  limit  computing  time  an 
artificial  interface  was  introduced  between  coarse  zoned  copper 
and  fine  zoned  copper.  The  result  showed  an  attenuation  of  150. 
kilobars  in  a  400  kilobar  shock  going  from  the  coarse  zoned 
copper  to  the  fine  zoned.  In  order  to  explain  this  spurious 
effect,  consider  the  neglected  terms  of  equation  1. 


n=a  2n-l  ni  '^n'j 


(A2mf  +  (_  Aim  f 


Am  +  A  m 
1  2 


The  largest  of  these  terms  (assuming  the  series  to  converge 
sufficiently  rapidly)  is  the  first: 


2  2! 


-  (Apm)2  +  (Ai.m)‘ 


A,m  +  A  m 
x  2 


E  =  -  (Ajin  -  A2m) 

4  dm2 
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Note  carefully  that  when  Ajffi  and  A  m  are  equal,  this  expression 
vanishes  and  the  error  is  reduced  to  third  order  rather  than 
second.  It  is  anticipated  that  the  third  order  terms  are  very 
small.  We  investigate  the  lead  term.  In  order  to  do  this,  we 
need  an  estimate  on  the  derivatives.  This  is  extremely 
difficult  to  obtain,  but  since  (5p/5m)  is  continuous  in  the 
smoothed  out  shock  front  (Richtmeyer  p203),  we  will  attempt  to 
estimate  dp/cta  and  32p/5m2  from  the  data  available  although  any 
estimates  of  higher  derivatives  probably  can  not  be  obtained  by 
the  method  we  will  employ.  The  pressure  is  initially  zero  and 
rises  to  .4  megabars  in  the  space  of  four  zones  -  2  to  the  left 
of  the  interface  and  2  to  the  right.  The  mass  increment  to 
the  left  of  the  interface  is  A^m  =  .0016  and  to  the  right  is 

.0001  (Agin) . 


To  estimate  the  second  derivative,  consider  that  5p/5m  is 
initially  zero  and  rises  to  the  average  value  listed  above 
in  a  space  of  one  zone  or  so  -  being  smaller  than  the  above 
value  in  the  first  zone  and  larger  in  the  second  so  that  the 
average  is  the  quantity  listed.  In  the  third  zone  the 
derivative  begins  to  decline  reaching  the  average  value  at 
the  end  of  the  third  zone  and  falling  off  to  zero  by  the  end 
of  the  fourth  zone. 

52p  op/ 5m  -  0 

5m2  2A2m 

5gp  0  -  5p/5m 
5m2  2  Aim 

Averaging  the  absolute  values,  we  ,et  a  crude  estimate  on  the 
maximum  value  of  52p/cim2: 
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d2p  ^  1  by/bm  by/bm 
5m2  2  2  A-jffl  2  Agin 

47  x  104 

5m2 


E  =  i-  i  £fE  (A.m  -  A.m) 

Si  2  5  m2  1  2 


=  i  (.47  x  104)  (  .0015) 

4' 

*•  180 

which  is  on  the  same  order  as  the  term  by/bm  which  is  being 
calculated.  Hence  near  a  shock  front,  a  mass  mismatch  can 
introduce  a  serious  error.  The  calculations  listed  are  by- 
no  means  sound,  they  are  merely  an  estimate  of  possible 
error.  At  the  present  writing,  an  attempt  is  being  made  to 
eliminate  the  difficulty,  but  results  are  not  yet  available. 

Another  question  that  frequently  arises  concerns 
equations  of  state.  As  is  pointed  out  in  the  text,  the 
confusion  arises  because  the  Hugoniot  is  widely  used  as 
an  equation  of  state  but  it  is  really  an  equation  of  process. 
However,  as  pointed  out  in  section  IV  (P2.6),  the  Hugoniot  is 
a  useful  approximation. 
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APPENDIX  A 

DERIVATIVE  OP  THE  JACOBIAN  DETERMINANT 


DERIVATIVE  OF  THE  JACOBIAN  DETERMINANT 


Let 


^  eijk  eenn  ^ie  ^km 
o# 

where  Tie  =  ~i.  and  the  are  the  familiar  alternating 
symbols  of  tensor  analysis  ( see  note  at  end) . 

SJ  _  1  STie  m  m 

ST'rS  3i  £iik  6emn  Tjm  Tkn 

-i  St, 

+  ~  eiik  eemn  r~^  Tle  Tkn 

3* 


i-  ST 

3i  eiJk  £emn 


kn 


rs 


T  T 

le 


3,  ^eUk  €eran  5ir  5es  Tjm  Tkn 

+  £ijk  eemn  8jr  5ms  Tie  Tkn 
+  6ijk  eemn  8kr  8ns  ^ie 


3!  ^€rJk  €smn  TJm  Tkn  +  6irk  €esn  Tie  Tkn 


+  eijr  Sens  Aie  AjmJ 


3]  Ssran  Tjm  Tkn  +  erik  esen  Tie  Tkn 

+  e  s  'T1  T  1 

rij  sem  ie 


Renumbering  the  dummy  indices  of  summation  we  obtain: 


So  that 


T  .  Fi  sj 

J.ps  -  - - 

LJ  5TrsJ 


=  5, 


pr 


or 


i§L-.  ppsl  -1 

j  L  J 


daw  ctoc  cta^, 

But  |  T  |  is  obviously  — £  since  — £  — £  =  5 


~  "l-i 

T 

rs 


Sx. 


eta  bxn 
r  p 


sp 


Hence,  as  claimed  in  the  body  of  the  report: 


■  -  =  t  ~1 

KI  r‘ 


J 


and 


T 


rs 


-l  ,  H 
Sx„ 


NOTE:  The  summation  convention  is  employed  throughout  so  that 
aij  "b jk  is  summed  from  J  =  1  to  3.  The  Kronecker  5,  5ltj,  has 
the  value  1  if  i  =  j  and  0  otherwise.  Obviously  bjj  a^  =  a^^ 
for  any  tensor  a.  Further  =  5u  +  622  +  533  =  3.  The 
alternating  symbol  e1Jk  is  zero  if  any  2  indices  (i, j,k)  are 
the  same.  It  is  one  if  i  j  k  is  an  even  permutation  of  1,2,3 
and  -1  for  an  odd  permutation. 
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